// Sparse neural network using multiple width nnSize neural layers with
// forward connected weight switching as the non-linearity. Uses 
// the fast Walsh Hadamard transform (WHT) as a connectionist device 
// (synthetic fully connected layer.)
// Uses a sub-random pattern of sign flips to prevent the first
// WHT producing a biased spectral response to natural data.
// Uses the P5 (Processing) JavaScript library.
// Sean O'Connor 20 September 2021 
// Public Domain
class SWNet {
  constructor(vecLen, nnSize, depth) {
    this.vecLen = vecLen;
    this.nnSize = nnSize;
    this.depth = depth;
    this.scale = 1.0 / sqrt(vecLen);
    this.params = new Float32Array(2 * nnSize * vecLen * depth);
    this.work = new Float32Array(nnSize);
    this.flips = new Float32Array(vecLen);
    for (let i = 0; i < vecLen; i++) {
      this.flips[i] = sin(i * 1.895) < 0 ? -this.scale : this.scale;
    }
    for (
      let i = 0, j = 0;
      i < this.params.length;
      i += 2 * nnSize, j = (j + 1) & (nnSize - 1)
    ) {
      this.params[i + j] = this.scale;
      this.params[i + j + nnSize] = this.scale;
    }
  }

  recall(result, input) {
    for (let i = 0; i < this.vecLen; i++) {
      result[i] = input[i] * this.flips[i];
    }
    let paIdx = 0; // parameter index
    for (let i = 0; i < this.depth; i++) {
      wht(result);
      for (let j = 0; j < this.vecLen; j += this.nnSize) {
        for (let k = 0; k < this.nnSize; k++) {
          let x = result[j + k];
          let idx = x < 0 ? paIdx : paIdx + this.nnSize;
          for (let m = 0; m < this.nnSize; m++) {
            this.work[m] += x * this.params[idx + m];
          }
          paIdx += 2 * this.nnSize;
        }
        for (let k = 0; k < this.nnSize; k++) {
          result[j + k] = this.work[k];
          this.work[k] = 0;
        }
      }
    }
    wht(result);
  }
}

// Fast Walsh Hadamard Transform provide your own scaling
function wht(vec) {
  let n = vec.length;
  let hs = 1;
  while (hs < n) {
    let i = 0;
    while (i < n) {
      const j = i + hs;
      while (i < j) {
        var a = vec[i];
        var b = vec[i + hs];
        vec[i] = a + b;
        vec[i + hs] = a - b;
        i += 1;
      }
      i += hs;
    }
    hs += hs;
  }
}

// Sum of squared difference cost
function costL2(vec, tar) {
  var cost = 0;
  for (var i = 0; i < vec.length; i++) {
    var e = vec[i] - tar[i];
    cost += e * e;
  }
  return cost;
}

class Mutator {
  constructor(size, precision, limit) {
    this.previous = new Float32Array(size);
    this.pIdx = new Int32Array(size);
    this.precision = precision;
    this.limit = limit;
  }
  mutate(vec) {
    for (let i = 0; i < this.previous.length; i++) {
      let rpos = int(random(vec.length));
      let v = vec[rpos];
      this.pIdx[i] = rpos;
      this.previous[i] = v;
      let m = 2 * this.limit * exp(random(-this.precision, 0));
      if (random() < 0.5) m = -m;
      let vm = v + m;
      if (vm >= this.limit) vm = v;
      if (vm <= -this.limit) vm = v;
      vec[rpos] = vm;
    }
  }
  undo(vec) {
    for (let i = this.previous.length - 1; i >= 0; i--) {
      vec[this.pIdx[i]] = this.previous[i];
    }
  }
}

// Test with Lissajous curves
let c1;
let ex = [];
let work = new Float32Array(256);
let parentCost = Number.POSITIVE_INFINITY;
let parentNet;
let mut;

function setup() {
  createCanvas(400, 400);
  parentNet = new SWNet(256, 8, 2);
  mut = new Mutator(5, 25, parentNet.scale);
  c1 = color("gold");
  for (let i = 0; i < 8; i++) {
    ex[i] = new Float32Array(256);
  }
  for (let i = 0; i < 127; i++) {
    // Training data
    let t = (i * 2 * PI) / 127;
    ex[0][2 * i] = sin(t);
    ex[0][2 * i + 1] = sin(2 * t);
    ex[1][2 * i] = sin(2 * t);
    ex[1][2 * i + 1] = sin(t);
    ex[2][2 * i] = sin(2 * t);
    ex[2][2 * i + 1] = sin(3 * t);
    ex[3][2 * i] = sin(3 * t);
    ex[3][2 * i + 1] = sin(2 * t);
    ex[4][2 * i] = sin(3 * t);
    ex[4][2 * i + 1] = sin(4 * t);
    ex[5][2 * i] = sin(4 * t);
    ex[5][2 * i + 1] = sin(3 * t);
    ex[6][2 * i] = sin(2 * t);
    ex[6][2 * i + 1] = sin(5 * t);
    ex[7][2 * i] = sin(5 * t);
    ex[7][2 * i + 1] = sin(2 * t);
  }
  textSize(16);
}

function draw() {
  background(0);
  loadPixels();
  for (let i = 0; i < 100; i++) {
    mut.mutate(parentNet.params);
    let cost = 0;
    for (let j = 0; j < 8; j++) {
      parentNet.recall(work, ex[j]);
      cost += costL2(work, ex[j]); // autoassociate
    }
    if (cost < parentCost) {
      parentCost = cost;
    } else {
      mut.undo(parentNet.params);
    }
  }
  fill(c1);
  for (let i = 0; i < 8; i++) {
    for (let j = 0; j < 255; j += 2) {
      set(25 + i * 40 + 18 * ex[i][j], 44 + 18 * ex[i][j + 1]);
    }
  }
  for (let i = 0; i < 8; i++) {
    parentNet.recall(work, ex[i]);
    for (let j = 0; j < 255; j += 2) {
      set(25 + i * 40 + 18 * work[j], 104 + 18 * work[j + 1]);
    }
  }
  updatePixels();
  text("Training Data", 5, 20);
  text("Autoassociative recall", 5, 80);
  text("Iterations: " + frameCount, 5, 150);
  text("Cost: " + parentCost.toFixed(3), 5, 170);
}

